Reducing the Variance

Assaraf and Caffarel have published a couple of papers (JCP 119, 10536 and JCP 113, 4028) about reducing the variance for general estimators in QMC. The VMC energy has the well known zero variance property - the variance goes to zero as the trial wavefunction approaches an eigenstate. Through a "renormalization" , the authors show how to extend this property to quantities other than the energy.

The basic approach seems to be to add a term (with zero average) to the bare estimator. Then that added term can be optimized to reduce the variance. Their derivation starts with the Hellmann-Feynman theorem, but this seems confusing to me.

Why not start with the new estimator, O tilde = O + B. ("O tilde" and "O" are notations from the paper. "B" is my own notation for representing all the added terms) Then we demand < B > = 0, and use a Lagrange multiplier to enforce it while optimizing the variance (<(O+B)^2> - < O+B>^2)?

Their paper gives a particular form for B, that ensures the < B>=0 constraint is met (or at least that < B > vanishes as the trial wavefunction approaches the ground state).

Can this approach be used to reduce the variance in classical Monte Carlo as well? In principle, B could optimized in that case as well - the hard part is finding a special form as in Quantum Monte Carlo.

Added thoughts (12/14/2003): The primary application of this method is to forces, which are plagued by an infinite variance in QMC. This looks a like good way to improve computations of the force.

The paper looks at average values and variances as functions (or at least the polynomial order) of delta psi (the diference between the exact ground state and the trial wavefunction). For general (non-energy) estimators, the bias (or error) in the average value is linear in delta psi, and the variance is constant. Using their procedure, the error is quadratic and the variance is linear. It appears the added terms (B) cancel out the leading error terms in the estimator. The question is then, why stop at the leading order? Can a similar process by applied to cancel out the quadratic term in the average ??

Free Energy Chronicles

Free energy calculations are annoyingly difficult. Computing averages with Metropolis is relatively easy,
but the free energy is hard to get. The usual method is thermodynamic integration over a path from a state of known free energy to the desired state. (Maybe it's the ensemble averages that are anomalously easy, making other quantities seem difficult by comparison) In addition, some other aspects about free energy calculations are bothering me as well.

The typical example for explaining Monte Carlo integration is to compute the area of an irregularly shaped field or pond by throwing rocks at it. But then computing the area (volume, free energy, whatever) is actually quite difficult in practice. Of course in this example, one has to compare the rocks falling in desired area to those falling in an enclosing area that can be computed analytically, so it's still computing the difference between two states.

Another issue that bothers my intuition is that if I were to wander around in a field I would have some idea how big the field is based on my wanderings. In other words, why can't we get the area from properties of the random walk. This could simply be pushing the analogy too far. A person wandering is not undergoing a truly random walk. They would likely walk in single direction until hitting an edge and then change direction. So imagine riding an out-of-control Segway around a field - would you still get a sense of the size of the field.

Additionally, a person is able to see the edges. So imagine being in a corn field (I am in central Illinois, after all - others can imagine a forest, with very soft trees) on an out-of-control Segway. Would it still be possible to get an idea of how big the field is based on the path taken and how often you bump into the edge?

So I'd like to look into some of the properties of random walks - how much information can be extracted. Are there theoretical limits? Can we change the random walk or come up with different estimators that would let us compute the area?

In order to investigate concrete but simple systems, I'm thinking of a series of systems to make the problem increasingly realistic for physical systems. The simplest and easiest to think about the area of an irregular region. More realistic, the area of the region is much less that the area of the enclosing square, and the region is long and snaky. Then go to higher dimensions, and eventually to a system of hard spheres. To be more generally applicable, the next step would be a general potential (think of a hilly field and we want to know the volume of soil is in the field)

Interesting notes

While searching for Monte Carlo information on the web, I ran across Jonathon Goodman's Monte-Carlo Course Materials. One interesting part was Alan Sokal's introductory notes near the bottom (see also this CiteSeer link). I liked these notes because they covered a practical method for computing the autocorrelation time.

Optimization

Recently there were two papers in the Journal of Chemical Physics related to QMC optimization. The first, Geometry optimization in quantum Monte Carlo with solution mapping: Application to formaldehyde, is concerned with optimization of nuclear positions using only energy information, and not gradient (force) information, since it is noisier and harder to obtain. They use statistical experiment design techniques to pick coordinates at which to perform the energy computations, and fit that data to a quadratic form. The the geometry was then optimized using that fitted form.

The second paper is concerned with optimizing the one-body part of the wavefunction: Transcorrelated method for electronic systems coupled with variational Monte Carlo calculation. The authors use a transcorrelated Hamiltonian (some of the correlation has been moved into the Hamilitonian), and compute the SCF solution for the one-body orbitals. They use regular VMC to optimize the Jastrow factors. Note that the two calcuations are coupled, and to reach a solution the procedure is to alternately optimize the one-body and Jastrow factors.

Monte Carlo on the Graphics Card

A recent short paper on The Graphics Card as Streaming Computer (see related presentation slides) got me wondering if Monte Carlo could be programmed onto a graphics card.

The proposed Protein Explorer from Japan, using a GRAPE type architecture, follows a similar programming model.

One limitation is the precision of the calculations. Single precision at best, some cards only have 8-16 bits. The penalty method or related methods may be useful for correcting errors due to finite precision. (See also Ball, Fink, and Bowler)

Another limitation is the size and complexity of the code that can execute on the graphics card. This might be more useful for simple classical Monte Carlo simluations that QMC.

I suspect this falls into the category of "cool ideas that would be fun to play with, but spending the time researching better algorithms is probably more effective".

Of Wavefunctions and Wikis

and mailing lists (but that didn't alliterate nicely in the title).

There's a new QMC Mailing list in town. Be the first on your block to sign up.

Also there's a QMC Wiki (A wiki is an easily editable set of community web pages) (No link yet - you'll have to sign up for the mailing list and look at the Sept. archives to get the address)

Both of these are brought to you by Alan Aspuru-Guzik.

Linear charged particle simulations

Today's paper is by Joerg Rottler and A.C. Maggs, entitled "A Continuum,O(N) Monte-Carlo algorithm for charged particles". There are some earlier papers in the references where this method was developed for lattice systems.

As near as I can tell, the algorithm works by writing the total energy in terms of the electric field rather than the electric potential. The electric field is put on a grid. Particle moves then only require local updates to the electric field grid.

The electrostatic interaction is generated as a result of minimizing the energy functional (functional of the electric field). The algorithm they give is supposed to recover the proper interaction by a Monte Carlo means. This is the part I don't quite understand yet.

It seems like a good example of exploiting Monte Carlo to reduce the complexity of a problem by solving it on average.